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O : ABSTRACT 

o , 

I In this paper we determine the onset point of secular instabihty for the non- 

1—5 I axisymmetric bar mode in rigidly rotating equilibrium configurations in the Post- 

^ I Newtonian approximation, in order to apply it to neutron stars. The treatment 

I is based on a precedent Newtonian analytic energy variational method which we 

^ I have extended to the Post-Newtonian case. This method, based upon Landau's 

OO I theory of second-order phase transitions, provides the critical value of the ellip- 

^ I sold polar eccentricity e at the onset of instability, i.e., at the bifurcation point 

^ I from the axisymmetric Maclaurin to the triaxial Jacobi ellipsoids, and it is valid 

^ I for any equation of state. The extension of this method to Post-Newtonian fluid 

^ I configurations has been accomplished by combining two earlier orthogonal works, 

' specialized respectively to slow rotating configurations but with arbitrary density 

^ I profile and to constant mass density but arbitrarily fast rotating ellipsoids. We 

I also determine the explicit expressions for the density functionals which allow the 

^ I generalization of the physical quantities involved in our treatment from the con- 

[ stant mass density to an arbitrary density profile form. We find that, considering 
homogeneous ellipsoids, the value of the critical eccentricity increases as the stars 



^ \ become more relativistic, in qualitatively agreement with previous investigations 

but with a less marked amount of such an increase. Then we have studied the de- 
pendence of this critical value on the configuration equation of state. Considering 
polytropic matter distributions, we find that the increase in the eccentricity at 
the onset of instability with the star compactness is confirmed for softer equations 
of state (with respect to the incompressible case). The amount of this stabilizing 
effect is nearly independent of the polytropic index. 

Subject headings: gravitation — instabilities — relativity — stars: interiors — 
stars: neutron — stars: rotation 
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1. Introduction 

The determination of the onset point of instabihty for nonaxisymmetric modes in rapidly 
rotating equihbrium configurations is a classic problem. In particular the m = 2, so-called 
"bar mode" has long been studied due to its relationship with the dissipation mechanisms 
of viscosity and gravitational radiation. In fact, there is a large number of astrophysical 
situations in which this instability may appear: the coalescence of a binary neutron star in 
a single, rapidly rotating object (Baumgarte et al. 1998); the core collapse in a massive, 
evolved star or the accretion-induced collapse of a white dwarf (Lai & Shapiro 1995); the 
coalescence of a white dwarf binary into a progenitor of Type la supernovae (Iben & Tutukov 
1984; Yungelson et al. 1994) or of isolated millisecond pulsars (Chen & Leonard 1993); the 
accretion and spin-up of a neutron star (NS) in an X— ray binary system (Chandrasekhar 
1970; Friedman & Schutz 1978); the implosion to a black hole of a supramassive neutron 
star (SMNS) after the spin-up phase (Salgado et al. 1994); the nonexplosive core contraction 
of a rapidly rotating massive star (a "fizzler", Hayashi, Eriguchi & Hashimoto 1999). And 
not only situations of coUisional systems: the bar formation in self-gravitating coUisionless 
galaxies in purely rotational equilibrium can also be studied as an application of this classic 
problem (Hohl 1971; Ostriker & Peebles 1973; Binney & Tremaine 1987). 

In Newtonian theory a considerable amount of work has been done, and a number 
of results are well established. Chandrasekhar (1969a), using the tensor virial formalism, 
found an exact analytic solution for the equilibrium shape and stability of an incompressible, 
homogeneous, rigidly rotating fluid configuration. In this case, the equilibrium shape in 
axisymmetry is a Maclaurin spheroid. Nonaxisymmetric instabilities develop in spinning 
spheroids when the ratio K/ | 14^ | of the rotational kinetic to the gravitational potential 
energy becomes sufficiently large. At the critical value K/ \ W 0.1375 the equilibrium 
sequence of Maclaurin configurations bifurcates in two different branches of triaxial equilibria, 
the Jacobi and Dedekind ellipsoids. Since the Maclaurin spheroids are dynamically unstable 
only for K/ \ W \> 0.2738, the bifurcation point is dynamically stable. However, in the 
presence of a dissipative mechanism such as viscosity or gravitational radiation, this point 
becomes secularly unstable to the 171 — 2 bar mode. 

After these findings, a number of efforts have been devoted to their extension to more 
reahstic, compressible fluids. Modeling the fluid with a polytropic equation of state (EOS), 
it was found that for strictly rigid rotation bifurcation to triaxial conflgurations can exist 
only when the polytropic index is less than the critical value n = 0.808 (Jeans 1919. 1928; 
James 1964; Tassoul & Ostriker 1970). This is because of the dynamical constraint that 
the angular velocity at the bifurcation point must be lower than the limiting value at which 
the centrifugal force balances the gravitational force at the equator ("mass-shedding" limit). 
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and this happens for equations of state sufficiently stiff. Later Ipser & Managan (1985) 
demonstrated that, while in general the m = 2 Jacobi and Dedekind bifurcation points do 
not have the same location along axisymmetric rotating sequences, when considering poly- 
tropic axisymmctric sequences with uniform rotation it is found that these two bifurcation 
points have indeed the same location, as in the incompressible case. Lai, Rasio & Shapiro 
(1993) obtained the same result constructing triaxial ellipsoid models of rotating polytropic 
stars in Newtonian gravity and using then an ellipsoidal energy variational method. This 
approach was originally introduced by Zeldovich & Novikov (1971) for the axisymmetric 
case and is also illustrated in Shapiro & Teukolsky (1983). Generahzing this approach to 
the triaxial case, Lai, Rasio & Shapiro (1993) were able to construct equilibrium sequences 
for compressible analogs of most classical incompressible objects, such as the Maclaurin, 
Jacobi and Dedekind ellipsoids. From these equilibrium and dynamical ellipsoid models, Lai 
& Shapiro (1995) confirmed that also in the compressible case the Maclaurin configurations 
bifurcate at K/ \ W |~ 0.1375, independently of the polytropic index n, and that a dy- 
namical bar mode instability sets in at K/ \ W 0.2738. In this paper it is also pointed 
out that the limit n — 0.808 imposed by mass-shedding in strictly uniformly rotating stars 
can be by-passed with a slight amount of differential rotation, which can in principle inhibit 
mass-shedding without changing the global structure of the star. On the other hand, a num- 
ber of authors have constructed numerical models of rotating equilibrium stars, using both 
polytropic (Bodenheimer & Ostrikcr 1973; Ostriker & Bodenheimer 1973; Managan 1985; 
Imamura, Friedman & Durisen 1985; Ipser & Lindblom 1990) and more realistic equations of 
state for white dwarfs and NSs {e.g., Ostriker & Tassoul 1969; Durisen 1975; Hachisu 1986; 
Bonazzola, Frieben & Gourgoulhon 1996). 

An analytic result independent of the EOS has been obtained by Bertin & Radicati 
(1976; hereafter BR). Using Landau's theory of second-order phase transitions (Landau & 
Lifshitz 1967) they found that the transition from the axisymmetric Maclaurin to the triaxial 
Jacobi sequence always corresponds to K/ \ W \^ 0.1375, if one assumes that the density 
is constant over ellipsoids with constant eccentricities and that both the internal energy 
and the enthalpy are independent of the shape of the rotating fluid. A formal treatment 
of the correspondence between the second order phase transition and the Maclaurin-Jacobi 
bifurcation is also presented in Hachisu & Eriguchi (1983). 

All these results may be characteristic of Newtonian physics and an gravitational 
force. For NSs and relativistic objects, the analysis of equilibrium and stability must be 
based necessarily on general relativistic models. If the structure of rotating axisymmetric 

stars in general relativity has been investigated numerically by a number of authors {e.g., 
Cook, Shapiro & Teukolsky 1992, 1994a,b, 1996; Salgado et al. 1994), the first fully rel- 
ativistic perturbation computations of nonaxisymmetric instabilities have been published 
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only recently by Stergioulas (1996) and Stergioulas & Friedman (1998), who focused on the 
instabihties driven by gravitational radiation. An earlier numerical investigation of the ef- 
fects of relativity on the viscosity- driven bar mode instability was carried out by Bonazzola, 
Frieben & Gourgoulhon (1996), and these results have been corroborated by a more detailed 
analysis (Bonazzola, Frieben & Gourgoulhon 1998). Stergioulas & Friedman (1998) find that 
relativistic models are unstable to nonaxisymmetric modes for significantly smaller degrees 
of rotation than for corresponding Newtonian models. The destabilizing effect of relativity 
is most striking in the case of the m — 2 bar mode, which can become unstable even for soft 
polytropes of index n < 1.3, to be compared with the Newtonian critical value n — 0.808 
reported before. This behaviour is in agreement with the result of a numerical investigation 
made by Yoshida & Eriguchi (1997) in the "relativistic Cowling approximation", in which 
all metric perturbations are omitted. However, the results of Bonazzola, Frieben & Gour- 
goulhon (1996), concerning the effects of general relativity on the viscosity-driven bar mode 
instability, seem to suggest the opposite effect. The critical polytropic index for the bar mode 
instability becomes lower, but only slightly, than the Newtonian value as the configuration 
becomes increasingly relativistic, reaching a value n ~0.71 for very relativistic objects. This 
behavior suggests that relativistic effects tend to stabihze the configurations. In noting this, 
Stergioulas & Friedman (1998) concluded that, in general relativity, the onset point of the 
viscosity- driven and the gravitational radiation-driven m = 2 modes may no longer coincide 
as they do in Newtonian theory, and that the effect of relativity seems to be very different 
in the two cases. Yoshida ct al. (2002) have recently found, again in the Cowling approxi- 
mation, that a further destabilizing effect in relativistic configurations is due to differential 
rotation. 

But from the analytical point of view, the general relativistic treatment of a nonax- 
isymmetric instability is difficult. Apart from the problem of solving Einstein's equations 
without any presupposed symmetry, in the relativistic regime emission of gravitational waves 
must be taken into account. However, a first step towards the evaluation of general relativis- 
tic effects can be made by examining the problem in the so-called "Post-Newtonian (FN) 
approximation" , where gravitational radiation can be neglected. To obtain this level of ap- 
proximation, the metric and the stress tensor are expanded as sums of terms of successively 
higher order in the expansion parameter c"^, where c is the velocity of light, while each 
equation is decomposed into a series of equations of successively higher order in c~^. The 
first FN correction will refer to the terms that are (9(c^^) (z.e., 0{GM/{c^R)), where M is 
the configuration mass, R a length scale of the problem and G the gravitational constant) 
higher than the corresponding Newtonian terms in this expansion. Gravitational radiation 
enters only at 2.5 FN level and higher. 

FN effects on the equilibrium of uniformly rotating, homogeneous objects have been 
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investigated by Chandrasekhar (1965a,b, 1967a,b,c, 1969b) using the tensor virial formal- 
ism, and even the Post-Post-Newtonian (PPN) corrections have been considered with this 
method (Chandrasekhar & Nutku 1969). In these works integral expressions for global con- 
served quantities are obtained, but no explicit formulae are given for the PN corrections to 
the rest mass, angular momentum and binding energy. A simpler formalism has been pro- 
posed only recently by Shapiro & Zanc (1998; hereafter SZ), who extended the Newtonian 
treatment of Lai, Rasio & Shapiro (1993) to PN gravitation. Considering incompressible, 
rigidly rotating bodies, they were able to construct equilibrium sequences of constant rest 
mass deriving the analytic functionals for the main global parameters characterizing a ro- 
tating configuration, and provided for the first time analytic investigations of the location 
of the secular instability point in general relativity. Their result is that the value of the 
ratio K/ \ W \^ defined invariantly, at the onset of bar mode instability increases as the 
stars become more relativistic, z.e., increases with the compactness parameter GM / {c^R)^ 
being K/ \ W |~ 0.1375 only in the Newtonian limit GM/{c^R) = 0. Since higher degrees 
of rotation are required to trigger a viscosity- driven bar mode instability as the stars be- 
come more compact, the effect of general relativity is to weaken the instability, at least to 
PN order. This behavior, consistent with Bonazzola, Prieben & Gourgoulhon (1996), but 
contrasting that found by Stergioulas & Friedman (1998), supported the suggestion that in 
general relativity nonaxisymmetric modes driven unstable by viscosity no longer coincide 
with those driven unstable by gravitational radiation. 

In this paper we wish to investigate analytically the location of bar mode instability 
points in rotating equilibrium PN configurations for arbitrary equations of state. The paper 

is organized as follows. We begin with a brief review of the second-order phase transition 
method of BR in §2, and then in §3 we start extending this method to PN configurations, 
obtaining the expression for the PN total energy (cqs. (54)-(58)). In §4 we determine the 
general expressions of the density functionals necessary to model any EOS (cqs. (62), (69), 
(77), (78), (82), (86), (91), (93), (96)). In §5 we give the complete treatment for the analytic 
determination of the onset point of bar mode instability, and in §6 we finally evaluate this 
point for various equations of state. In §7 we report a discussion of our findings and finally 
in §8 the conclusions of our work. 



2. The Newtonian treatment of Berlin &; Radicati 

As previously reported, the treatment made by BR of the nonaxisymmetric instability 
which leads from the axisymmetric Maclaurin to the triaxial Jacobi sequence is based on 
Landau's theory of second-order phase transitions (Landau & Lifshitz 1967). Second-order 
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phase transitions occur in crystals when, as the temperature decreases, the invariance group 
of the crystal suddenly reduces to one of its subgroups. In the phase with lower symmetry 
a new observable, the "order parameter" ^, which vanishes in the symmetrical phase, is 
necessary to describe the state of the system together with the thermodynamical variables 
such as the pressure P and the temperature T. In general a second-order phase transition 
occurs along a line in the (P, T) plane which divides it in two regions corresponding to 
different symmetries: in region I and on the transition curve the order parameter ^ vanishes, 
and we have higher symmetry, while in region II ^ > and there is lower symmetry. 

Consider now the total energy E as a function of the thermodynamical variables: entropy 
S, volume V, angular momentum J, and of the order parameter ^. Expanding E in powers 
of ^ in the neighborhood of ^ = 0, one gets: 

E = Eo{S, V, J) + iE,{S, V, J) + eE2{S, V, J) + eEs{S, V, J) + ... (1) 

For the equilibrium of a physical configuration with such a total energy, it must be 
dE/d^ — 0. Now there are two possible solutions: one with ^ = and the other when ^ > 0. 
In this latter case, in which obviously E\ — 0, if it is also £'3 = the change in stability is 
defined by the condition: 

E2{S,V,J)=^ (2) 

To discuss the symmetry breaking in the case of a self-gravitating fluid, BR considered 
that in the symmetrical phase the shape of the fluid is an ellipsoid with polar eccentricity 
< e < 1 and that, as a result of the symmetry breaking induced by the instability, the 
system acquires an equatorial eccentricity ^ > 0. In this notation the Maclaurin sequence 
is thus characterized by ^ = 0, while along the Jacobi sequence we find the solutions with 

e>o. 

BR make the following general assumptions: 

(i) the internal energy and the enthalpy are independent of the shape of the rotating fluid; 

(ii) the density is constant over ellipsoids with constant eccentricities ("ellipsoidal approxi- 
mation" ) . 

Both these assumptions are discussed in Zeldovich & Novikov (1971), and BR show that 
as a consequence the critical value of the ratio Kj | 1^ |, or equivalently of the eccentricity e, 
where the Jacobi sequence branches off from the axisymmetric Maclaurin sequence, is given 
by the transition of E2 from positive to negative sign, ie., by the condition (2). This result 
is not restricted to a particular EOS. 

In order to determine this critical value, BR write the total energy as the sum of the 
gravitational, rotational and internal energies: E — W + K + U. Assumptions (i) and (ii) 
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imply that the internal energy U is independent of the shape while the gravitational and 
rotational energies are ^: 

W = --(-j ^/3[p]^(e,0 (3) 

thus depending on the fluid shape. Here M is the mass and 7[p], are functionals of the 
density p which reduce to 1 for constant p. The functions / and g are: 

. ^'^---)':°(^-^)'^' (5) 



Minimizing the total energy E with respect to the volume V (scalar virial equation) 
they obtain expressions for the first few terms in expansion (1) for the total energy in the 
neighborhood of ^ = 0, as a function of the eccentricities e, S,- By factoring out the part of 
W independent of e, ^ (and calling it Wq) and using the notation = df /dx, the result is: 

E 1 

Tif = -7{g^^ + '^^ge^ + ^^gee + Qe^^ + ^ge^e) (8) 
Wo 4 

where: 

\9Us,v,j 

(the variables that appear in the subscript of this latter definition are considered constant) . 
All the derivatives must be calculated at ^ = 0. 

Requiring then the validity of the two equilibrium conditions for the total energy of the 
ellipsoidal configurations, dE/de = and dE/d^ = 0ate7^0,^ = 0, BR first verify that 
£■1 = and then calculate the value of e for which the condition £'2 = is satisfied. This 
value (ec = 0.81267) corresponds to a critical ratio K/ \W \— 0.13752. 



^In BR's original expressions for W and K the numerical factors were omitted. Here we restore them. 
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3. The total energy in the PN approximation 

One of the goals of this paper is to push to PN order BR's version of Zeldovich and 
Novikov's (1971) energy variational method, which we briefly described in the previous sec- 
tion. Like Zeldovich and Novikov, and every author since then, we shall make the assumption 
that the density is constant on ellipsoids of fixed eccentricities e, { (we follow BR's choice of 
symbols) . We thus need to find exphcit expressions for the PN corrections to both the kinetic 
and the potential energy of a fluid conflguration, and for arbitrary eccentricities. These ex- 
pressions do not exist in the literature: PN corrections have been given by Bisnovaty-Kogan 
and Ruzmaikin (1973; hereafter BKR) and SZ, but while SZ have considered arbitrary eccen- 
tricities, but uniform density, BKR considered arbitrary density proflles, but conflgurations 
with a small deviation from the spherical shape. 

To obtain the explicit expressions for such PN corrections we will combine both BKR's 
and SZ's works. Therefore it is useful flrst to summarize the results of these two papers. 

3.1. Bisnovaty-Kogan &; Ruzmaikin's work 

In their paper, BKR investigate the stability of rotating supermassive stars (SMS), 
i.e., those with M > 10^ M©, by adding to the usual expressions for the full mass-energy, 

rest mass and angular momentum the deviations arising from the flrst and second order 
PN corrections of general relativity in stationary rotating conflgurations. However, in this 
particular approximate energy variational method, slowly rotating (Hartle 1967) SMS are 
considered, and thus the kinetic energy terms due to rotation appear as corrections of the 
same order of the flrst PN corrections to the gravitational and internal energies. Similarly, 
the first PN corrections to the rotational kinetic energy result of the same order of the PPN 
corrections to the gravitational and internal energies. 

BKR start by choosing a metric which allows them to integrate the field equations more 
easily, and to derive from these latter expressions for the total mass-energy and angular 
momentum of a stationary rotating configuration. In a spherical coordinate system R, 9, 
the element of four-space of this metric takes the form: 

ds'^ = e^icdt - gR^ sin^ ed(j)f - e^{dR^ + R^dO"^ + e^R^ sin^ ed(i)'^) (10) 

where the independent functions A, /x, g are chosen to satisfy Einstein equations. When 
g. fj, — > 0, this metric reduces to the spherically symmetric isotropic one (Landau & Lifshitz 
1971). In this case R corresponds to the so-called "isotropic coordinate". 

Moving to the next point, BKR expand the chosen metric with a power series in c~^, and 
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calculate the total mass-energy and angular momentum in the PPN approximation. Then, 
in order to calculate the first and second order corrections to the Newtonian total energy, 
they move from the coordinate R to the Newtonian radius r. For this transformation, BKR 
use the relationship: 

where di and d2 are functions of the mass m and radius r. In their work, BKR report the 
following expressions of di and c?2 correct to the PPN order: 
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where ro is the Newtonian radius of the sphere, u is the internal energy per unit mass and 
Q the angular velocity. 

A by-product resulting from the equations given in BKR, which we will use in the 
following of the paper, is the expression of the rest mass of the sphere Mq in terms of 
coordinate R. At PN order, this is: 
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After the transformation of coordinates, BKR calculate the PN and PPN corrections 
El and En to the Newtonian energy E^- The results are the following ^: 




^In BKR's original expression for En we have found a few typos. We report here the corrected form. 
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where P is the pressure and p the density profile of the rest mass. In these expressions, 
integration is carried out over the whole mass of the star, and the limits of the interior 
integration go from the centre to the actual m or r. 

BKR consider also a correction Bob caused by the stellar oblateness: 

2 vn / pro jy^d^ 
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where the value of a defines the degree of oblateness. About this value, BKR just note that 
it is of the same order of the ratio 2GM/ {(?R) between the Schwarzschild radius and the 
physical radius R of the star. However, since a gives a measure of the stellar oblateness, and 
the oblateness is caused by the rotation of the star, such a parameter must be also related 
to VL. Therefore, in the case of slow rotation treated by BKR, the correction Eoh is of the 
PPN order. 

Finally, BKR find the PN-corrected expression of the angular momentum J specialized 
to the case of constant f2: 
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3.2. Shapiro Zane's work 

As we already briefly mentioned in the Introduction, SZ construct analytic models of 
incompressible, uniformly rotating stars in PN gravity in order to evaluate their stability 
against nonaxisymmetric bar modes. For this, an energy variational principle is employed, 
its equations being exact at PN order. Contrary to BKR's work, this analysis is not restricted 
to slow rotation, whereby one requires {R^/{GM)y^^ Q <S 1, but arbitrarily fast rotation 
is allowed, so that Jl^ is permitted to reach ~ (GM/R^) and stars can suffer considerable 
rotational distortion. However, this particTilar energy variational method is valid only for 
constant density po- This latter condition implies that the internal energy vanishes and the 
Newtonian total energy is just given by: 

E^W + K (19) 

In the choice of the metric, SZ adopt a 3 + 1 ADM sphtted form (Arnowitt, Deser & 
Misner 1962) to solve Einstein's equations of general relativity. A subset of these equations 
reveals to be well suited to numerical integration in the case of strong-field, three-dimensional 
configurations in quasi-equilibrium. Moreover, the adopted equations are exact at PN or- 
der, where they admit an analytic solution for homogeneous ellipsoids. The most general 
expression for this kind of metric is: 

ds'^ = -a^dt^ + -fij (dx' + (3'dt) (dx^ + (3Ht) (20) 

where a and /5* are the lapse and shift functions, respectively. Then SZ choose a "conformally 
flat" decomposition of the spatial metric: 

lij = (21) 

where ^ is the "conformal factor" and fij is the Euchdean metric in the adopted coordinate 
system. SZ use Cartesian coordinates Xi (i=l,2,3). 

At this point, SZ expand the metric and the stress tensor in terms of c~^, and decompose 
each ADM equation into a series of equations of different order in c~^. To work in the PN 
approximation, they retain only the Newtonian terms and those that are C(c~^) higher. 
After such an expansion, they evaluate the conserved quantities of total mass-energy M, 
total rest mass Mq and angular momentum J, flrst in the integral form and then performing 
the quadratures over the fluid volume, adopting constant density triaxial ellipsoids with 
semiaxes of the outer surface specified by the values aj(i = 1, 2, 3). 

Since also the energy of the fluid E — {M — Mq)c^ is a conserved quantity, they exphcitly 
report it in the integrated form. For our later applications, we write down here their resulting 
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expressions for E and J, reintroducing the gravitational constant G and the velocity of hght 
c (SZ use geometrized units with c — G — 1 throughout their paper): 

7 « fiM.ii=i-(l + 5^P3/>) (23) 
where /, h, gi, Pi are functions of the ellipsoid axial ratios: 

2/3 / N 2/3 



Ai 



(3 ^ ^-S) 



while is the angular velocity of the fluid system and Mc is the quantity: 

function of the "conformal radial coordinate" R. From this latter definition it is possible 
to note that R represents the radius of the spherical configuration with the same volume as 
the rotating one, and thus it can be considered as a "mean radius" . The PN relationship 
between the coordinate quantity Mc and the total baryon rest mass Mq is also given: 

The structure of expression (22) is particularly convenient for performing the required 
energy variational method, since the full dependence on the two axial ratios is contained in 
/, /i, and, for the PN contributions, in Qi and pi. They have also checked that their result 
agrees with the PN correction to Newtonian energy obtained by Shapiro & Teukolsky (1983) 
for nonrotating, homogeneous spheres (see §B2 in Appendix B to their paper). For these 
latter configurations, R corresponds to the conformal isotropic coordinate. 

In order to construct sequences of axisymmetric equilibrium models, SZ first rewrite the 
energy £^ as a function of J, using the relationship: 

^2„2 25/^2 / 5GMe 25/^2 / GMc ,\ 

^'^'^77^-^ l + oT^PsM ^T777?7-^ 1-5-^7^ M (27) 



MIR^ 4 V 2 d^R J M'^R^ A \ d^R ^ 

which derives from expression (23), and combining it with expression (22). At PN order, 
they obtain: 

where P123 = P12 — Pz- Then the equilibrium sequence is determined by minimizing E with 
respect to Ai and A2 holding constant Mq and J. 
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3.3. Energy functional for eirbitrciry configurations 

Now we will combine the orthogonal sets of results presented in the two previous sub- 
sections, so as to obtain the general expressions of PN corrections to both the kinetic and 
the potential energy of a fluid configuration with arbitrary density profile and arbitrary 
eccentricities e, We do it as follows. 

Considering the kinetic energy, we write it as: 

Ktot + K^r (29) 

where K^orr is the correction to the Newtonian kinetic energy, complete to all orders, not 
just the PN one. For obvious dimensional reasons, we must have: 

where J, M, V are the model angular momentum, mass and volume, respectively, and L[p\ is 
an adimensional functional {i.e., an application which associates a number to every function) 
of the density profile p{r). We now rewrite L[p] as: 

L[p]=L[po]^^ (31) 

where po is the density of the constant mass density model with the same total mass, volume 
and shape as the stratified model. Obviously, again for dimensional reasons: 

= g y^^) (32) 

The argument of q is obviously the ratio of the Schwarschild radius to the physical radius. 

We now introduce the fact that we are only interested in PN corrections. Clearly, we need 
to introduce the hypothesis that, as GMj (c^y^/^) 0, both functions p(a;) and g(,T) are 
analytic at x = 0. In other words, we are assuming that it makes sense to expand general 
relativity terms in powers of GM/(c^l^^/^), which seems innocuous enough. In this way, we 
find: 

where l{e, ^) is a function of the model shape only. The first, explicit term is the PN correction 
to the kinetic energy, for arbitrary shape and density profile. However, by specializing to the 
case p — constant, we now see that l{e,^) must be exactly the function determined by SZ, 
just rewritten in terms of the polar and equatorial eccentricities e and ^ via the relationships: 

A,^(l-e^)V3 ; A,^fi-|j (34) 
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Also, when the model is spherical, L[p\/L[po\ must be the functional determined by BKR's 
work. In summary, we take for the PN correction to the kinetic energy: 

An entirely similar argument yields AW, the PN correction to the potential energy. In 
this case, we have: 

WtOT + Wcorr (36) 

where the correction Wcorr, complete to all orders, is: 

Wcorr^^p{L[p],e,0 (37) 
Expansion of the general relativity terms in powers of GM/ (c^V^^^) gives: 

GM^ GM L[p] ,^({ GMV\ 

where h{e,C,) is the eccentricity-transformed shape function gi2 which appears in SZ's PN 
correction to the potential energy. In summary, we can write: 



AW = ^ (M) 



{h{e,i))sz (39) 

BKR 



However, there is still an obstacle to the straightforward determination of the explicit 
expressions for the PN corrections AW and AK. In fact, from the previous discussion 
of BKR's and SZ's earlier works, it comes out that their PN corrections are expressed in 
different radial coordinates, respectively Newtonian radius in BKR and conformal radius 
Rc in SZ. Thus, in combining the results of these two papers, we must be careful with this 
difference. 

We start by rewriting BKR's PN expression (18) for the angular momentum in the case 
of rigid rotation as: 

= - ^I^MMkH - \Mml.iA (i - gH^^) (40) 

This form is similar to that adopted by BR for the expressions of the gravitational and 
rotational energies W and K (see eqs. (3)-(4)), and in a similar way (j[p] and r[p] are 
functional of the density profile which reduce to 1 for constant p. In §4, when we will 
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determine the explicit form of all the density functionals introduced in this section, the 
origin of the factors that multiply each of them will become clearer. 

Then we consider the resulting expressions (15)-(16) of BKR for the PN corrections to 
the Newtonian total energy, omitting those related to the internal energy which vanish in 
SZ's case of constant matter density distribution. Keeping in mind that for slow rotation 
the kinetic energy enters only as a correction, as already noted in §3.1, we take the terms up 
to the 0{c~^) order and write this total energy in the form, valid for uniform rotation: 

Ebkr = -^G-^Pip] + -Mn'R%^[p] - ^^^^ + T^^^'<« P 41 
5 Kn 5 70 C^Kj^ 175 c^itjv 

where we have introduced other functionals of the density p. In particular, this functional 
(3[p\ is exactly the same of that appearing in eq. (3). 

Now we exploit eq. (40) to find: 

^ 5 Jbkr ( _ 17 GMr[p]\ 

2MRla[p\ V 63c2i?jva[p]y ^ ' 

and therefore the PN-approximated relationship: 

o2 o2 ^ 25 Jl^^ ( 34 GM r[p] \ 

Introducing this latter relationship in eq. (41) we can rewrite the PN total energy in terms 
of angular momentum J instead of constant angular velocity VL. We obtain: 

^^^^ - "5^i?^^f^] + 4Mi4^-70^^[^] 
, 1GJ|^ 1 /85 7[p]r[p] 23 A 

Since it is possible to write M — [4:7r/3)poR%, indicating with po the mean mass density over 
the whole configuration, we can reduce the latter equation to an expression dependent only 
on the Newtonian radius i?^ and the angular momentum. We get: 

1GJ|^ 1 /85 7[p]rH 23 \ 

At this point, since we are operating in the slow rotation regime, we can use BKR's 
relationship (11) between the Newtonian and the conformal (quasi- isotropic) radial coordi- 
nates. After a proper treatment of eqs. (12) and (13) for the "transformation functions" di 
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and (^2, evaluated at the star boundary, we can write: 

where other two density functionals, fi[p] and z/[p], have been introduced. Exploiting again 
eq. (43), we find the same relationship as a function of the angular momentum: 

5^^[^] - Y2-^WW^^]) ^^^^ 

In order to move to SZ's radial coordinate i?c, we invert this latter relationship obtaining, 
at PN order: 

6 GM 5 Jsz ^[P] 



Substituting eq. (48) into eq. (45) , we obtain the total energy in terms of the coordinates 
adopted by SZ. At PN order: 

3 f47r\\ 2r.5^r 1 A QGM . . 5 u[p] 
Esz ~ -7 hr GplRl(3[p\ 1 + --p^p[p\ ' 



5 / 3 ^ Jlz 7[p] A 6 GM . . , 5 z/[p] ^ 



+ 4 V47r J po^^ ^'[p] V ^ 5 c2i?/f^^ ^ 12 c2M2i?2 a2[p] 1 ^^^^ 



70 i,T ) ^p^^^^^^ +14^^ U^i;^ ^ ^"^^ 



and after expansion in the PN correction terms, and since now M — [4'k/3)poRI (^^^ defini- 
tion (25)), we find: 

Finally, in order to generalize this expression of the total energy to arbitrary fast rotation, 
and thus arbitrary eccentricities e, ^, we introduce SZ's shape functions in such a way that 
expression (50) is just the limit for e, ^ — > 0. The result is: 

175 1 ^85M*M |,^HP1-15*]MP1 + ?I«[P1V(.0 (51) 



536 c2i?3 (72[p] V63 a[p] 2^'-'^'^'^' /la-ja-la-j -^^ 
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where we have rewritten SZ's functions / and h in terms of e, ^ and named them g{e, ^) and 
/(e,0 since they are exactly the same of those defined in BR and reported in eqs. (5) and 
(6). For what concerns the shape functions h{e,^) and l{e,^), whose origin we have already 
discussed above, we present here their full expressions: 



27 
'28 
99 
56 



-4? 



1 - e2)2/3(l _ ^)2/3 
l_^)l/3 



A 



(1 



,1/3 



(1 - e2)2/3 



AIA3 



2(1 _g2)2/3 

1 - e2)V3 



A. 



(l-e2)4/3' 



(1 - ^2/3 



- 0'/' 

+ A,A,{l-0'^'{l-e'y/' 



(52) 



39 3 (l-e2)2/3 ■ 
28^^40(1-0'/^ \ 



+ f 



J 



33 (1-0'/' 



{A,+A,) 



(53) 



280(l-e2)2/3' 

The dimensionless coefficients Ai are given in Chandrasekhar (1969a); for their calculation 
in terms of standard incomplete elliptic integrals involving only the eccentricities e and 
see Appendix A. The integral J is that reported in Appendix C of SZ. To transform it into 
a function of e, we used the relationships (34) between the axial ratios Ai, A2 and the two 
configuration eccentricities. 

The works of BKR and SZ have a common domain of vahdity, that of slowly rotating, 
constant density models, for which all functionals reduce to 1, and the shape functions are to 
be computed for e, ^ 0. The numerical factor for the PN correction to the kinetic energy 
in BKR (eq. (51)) differs from that of SZ (eq. (28)): to wit, -457/63 vs -67/7, respectively. 
The same disagreement is to be found in the PN corrections to both the angular momentum 
(from eq. (112) for BKR and eq. (23) for SZ we find 1342/1575 and 82/35 respectively). All 
other terms, instead, coincide. Faced with this dilemma, we remark that the derivation of 
the coefficients in SZ is buried in heavy computations, most of which are not reported, and 
which we could not thus check. The computations of BKR, instead, are fully detailed, and 
have allowed us an independent rederivation and a step-by-step comparison, from which we 
deduced that their work is surely correct. At the same time, the computations by SZ satisfy 
identically cqs. (106)-(108), which arc the obvious PN generalization of the equilibrium 
conditions in the Newtonian regime, for ,^ ^ 0. It thus appears that SZ's equations have 
the correct low triaxality limit, except for an overall factor. We have thus decided to adopt 
SZ's shape functions, but not the overall factor of the PN terms, which we take instead from 
BKR's work. This surely leads to the correct first order terms, in the limit ^ — > 0, as discussed 
above. We could not check independently the next term in ^ of SZ's shape functions: we 
simply assumed that they are correct, and we shall check this assumption by comparing our 
results with numerical investigations (see §7). Should any further modification be required, it 
is perfectly clear how one should proceed: in fact, all of the density functionals to be derived 
here are independent of the overall factors and of the shape functions, and the derivation 



-18- 



of the critical eccentricity to the bar mode instabihty also remains unaffected. At most, 
marginal numerical differences may result. 



3.4. The final expression for tiie PN total energy 

We can now write the explicit expression of the rotating configuration total energy E 
for arbitrary density profiles and eccentricities, which we will use in the PN extension of 
BR's treatment of nonaxisymmetric bar mode instability. 

The total energy of the fluid will be of the form: 

E^W + K + U + AW + AK (54) 

Focusing on each single term of this latter expression, we have that the form of the Newtonian 
gravitational energy W is exactly that of eq. (3), but with the dimensional quantities referred 
to conformal coordinates. We repeat it here for completeness: 

The Newtonian kinetic energy is given by: 

while the PN corrections to these two different kinds of energy contributions are: 

14 /47r\ 5 G'^M^ ( 1 \ 

= ^(yj ^(6/?WMp] + n'^[p])Me,0 (57) 

AX = 175 /^47r>^ 1 (%^iW\p\ 



536 V 3 y c^V (t2[p] V63 a{p\ 

5 23 \ 

—^my{p\ - 157[P]/^[P] + W J ^(e' (58) 

Because the internal energy U is independent of the rotating fluid shape, it is not necessary 
to write its explicit form. 

In order to extend BR's method to PN configurations, we still must determine the 
explicit expressions for the density functionals that we have introduced in this chapter. This 
will be the aim of next section. 
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4. The expressions of the density functionals 



When, in their work, BR write the gravitational and rotational energies in the forms 
(3)- (4), they do not give the explicit expression of the newly introduced density functionals 
f3[p\ and 7[p], and we have not found these expressions in the literature. We thus obtained 
them, together with all the other density functionals introduced in the previous section, by 
means of the following argument. 

The possibility of writing the gravitational and rotational energy terms as in eqs. (3) 
and (4) is due to a theorem of dimensional analysis (the so-called H-theorem; see, e.g., 
Barenblatt 1996). The density functionals appear separated from the shape functions f,g, 
their expressions are independent of e and ^, and therefore they can be calculated in the 
simpler spherical case. Moreover, from §3.3 it is evident that also for the other functionals 
the determination can be made in this simple case, and all the results found by BKR can be 
exploited. 

The main property of such density functionals is that they generalize the expression 
of a particular physical quantity from the constant density form to the arbitrary density 
distribution form, keeping fixed the values of the other physical parameters, and reducing 
to 1 for constant p. Their general form will thus be obtained by the ratio of the physical 
quantity with which they are related, written for an arbitrary density distribution, and the 
expression of the same quantity in the case of constant density. 

The determination of these expressions can be made using Newtonian coordinates, since 
the Newtonian contributions to the total energy take the same form in both coordinate 
systems, and moreover any PN correction to the density functionals of the PN energy terms 
would be of PPN order and therefore not interesting in our treatment. 

In the rest of this section we will consider each density functional previously introduced 
and determine its explicit form. 

4.1. The density functional for the Newtonian gravitational energy 

The functional is given by the ratio of the gravitational energy W for an arbitrary 

density distribution p(r) and that for a constant density po, at fixed values for M and 
V . Indicating with R the Newtonian radius of the spherical star and with m(r) the mass 
contained in a sphere of radius r: 




(59) 
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we thus have: 



W^- J ^^I^p{r)dV = -IGn^G rp{r)dr r'^p{r')dr' (60) 
For constant density: 

W ^ -l^^GAB!' ^ -IgMI (61) 

Note that this latter expression is exactly the factor which multiplies the density functional 
P[p] in the first term of the right-hand side of eq. (41). In general, the dimensional factor 
which appear before the density functional of a physical quantity is given by the expression 
of that particular physical quantity in the case of constant density matter distribution. This 
is a consequence of the definition of a density functional. 

The ratio between the last two expressions gives our functional: 



4.2. The density functional for the Newtonian kinetic energy 

In the case of the functional j[p], the ratio between the two rotational energies for 
different mass density distributions must be evaluated keeping fixed also the value of the 
angular momentum J. Calling Q the uniform angular velocity of the configuration with 
arbitrary density distribution and cu that of the constant density configuration, we have that 
the rotational energy can be written as: 

2 J 3 Jo 

which for constant density becomes: 

K = —TTu^poR^ = -Muj^R^ (64) 
15 5 

and therefore we obtain the functional expression: 

^i^^^^^'f^'"*^^* '''' 

But the condition of fixed angular momentum enables us to find the relationship between 
the two angular velocities ^2 and cu. Since the rotational energy can also be written as 
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K — (1/2)7Q^, where I is the momentum of inertia, from eq. (63) we obtain the integral 
expression for this latter physical quantity: 

8 

I — -n r^p{r)dr (66) 
3 Jo 

which for constant density gives: 

/ = ^^PoR' = \mR' (67) 

Thus the condition that the angular momentum J = IQ, must be the same in both configu- 
rations implies the relationship: 

n PoR^ 



^ 5 r'^p{r)dr 
Introducing this ratio in eq. (65) the functional j[p] results: 



(68) 



7[P] = rr/^ (69) 



4.3. The density functionals for the PN gravitational and kinetic corrections 



Considering now the functionals S[p] and a[p], to compute them we can recall, as ex- 
plained above, BKR's expressions (15)-(16) for the energy corrections up to the PPN order 
for an arbitrary density distribution. Retaining only the 0{c~^) terms, we can write the PN 
gravitational and kinetic corrections, in the case of constant angular velocity, in the form: 



AW = — - - 



+ 



AK 



G"^ (I f m?{r-)p{r)dV f p{r)dV f m{r)p{r)dV 



m{r)p{r)dV 



m{r)rdr 



(70) 



m{r)p{r)dV f p{r)dV 



+ 



u{r)p{r)dV 



IG 
-2 



4 1 
3^ 

p(r)dV 



(^J r^p{r)dV^ -^J P^I^ J r^p{r)dV + J m{r)rp{r)dV 



J'-^Jrnir)rdr-J^^^}^J 



r dr 



(71) 



+ 



1^1 
3 c2 



u(r) + 2— ) r^p(r)dV 
P. 
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Substituting eq. (59) in both these last expressions, we obtain: 



AW 



647r' 



1 

'2 Jo 



p{r)dr ( / r'^p{r')dr' 



+ 



r rp{r)dr f r' p{r')dr' f r"^p{r")dr" 
Jo Jo Jo 



p{r) 



dr I r'^p{r')dr' I r'dr' I r"^p{r")dr" 



,G 



AK = ^TT^f^^S 



pR rr pR pr 

' u{r)rp{r)dr I r'^p{r')dr' + / rp{r)dr / u{r')r'^p{r')dr' 
Jo Jo Jo 



(72) 



pR pr pR pr pr' 

+ / r^p{r)dr / r''^p{r')dr' - 2 / rp{r)dr / r'dr' / r"'^p{r")dr" 
Jo Jo Jo Jo Jo 

^dr f r'^p{r')dr' f r'^dr' 
^ Jo Jo 



(73) 



+^^' J r^p{r)dr J 



u{r) + 2— 1 sin^d^ 
P. 



which for constant density (and thus u = 0) give, after some calculations: 

32 , 3 Q2^3 



AW 
AK 



315 

368 2 C 2 2 n7 



70 C2i?2 

23 GM^ 



2 d2 



(74) 

1575" c^'^"" " 175 cm" " ^^^^ 
In the computation of AK (in particular for the last integration) we used the Newtonian 
result for the pressure P at constant density po (see, e.g., Chandrasekhar 1965a). which is 
also reported in eq. (55) of SZ's paper in terms of Cartesian coordinates. We rewrite it here 
in spherical coordinates: 



L = \gpo{R'' - r^) + ^a;V2 sin^ 9 

Po 3 2 



(76) 



The last term on the right-hand side of this relationship is negligible in the calculation of 
AK because of the slow rotation approximation that we are adopting in order to determine 
the density functionals. 

From the above equations, we thus obtain the following expressions for S[p] and a[p]: 



630 



p{r)dr I / r' p{r')dr' 



R 



rp{r)dr I r'p{r')dr' I r"^p{r")dr" 
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+ / ^dr I r'^p{r')dr' I r'dr' I r"^p{r")dr" 











^0 



(77) 



315 



a\p] = 



2'kGpIW 

7 



R pr pR 

u{r)r p{r)dr I r'^p{r')dr' + / rp{r)dr I u{r')r'^p{r')dr' 



r'^p{r)dr^ ~ ^ J ^P^^^^^ J ^"^ Pi^')^^' 

+^ rr^p{r)dr f r'"" p{r')dr' - & ^ rp{r)dr f r'dr' f r"^p{r")dr" (78) 
5 Jo Jo Jo Jo Jo 



4.4. The density functional for the Newtonian angulcir momentum 

The density functional a[p] has been introduced in eq. (40) in order to generahze the 
Newtonian constant density angular momentum to the case of arbitrary density distribution. 
Now, from eq. (18) we see that the Newtonian angular momentum in this latter case is: 

J=-n f r^p{r)dV = -nQ [ r^p{r)dr (79) 
3 Jo 3 Jo 

which for constant density becomes: 

J = —T^upoR^ = -MuR^ (80) 
15 5 

The ratio between eqs. (79) and (80) gives for a[p\ the result: 

Using relationship (68) between the two different angular velocities Q and uj we find that in 
the case of this density functional the simple identity: 

a[p\ = 1 (82) 

is valid for any density profile of the rotating configuration. 
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4.5. The density functional for the PN correction to the angulcir momentum 



BKR's result reported in eq. (18) gives us the PN correction to the angular momentum, 
which we will call A J, in the form: 



AJ 



2Q 
3c2 



m{r)rdr 



m{r)rp{r)dV 



AG p{r)dV 



Introducing the explicit eq. (59) for m(r), this gives 



(83) 



AJ 



r p{r)dr 



uir) + 



P 



sin 9d9 



64 



R 



r p{r)dr / r' p{r')dr' 



(84) 



2 f^' 
+ - / rp{r)dr I r'dr' I r"^p{r")dr" 
•J Jo Jo Jo 



Considering a constant density, and using relationship (76) for the ratio P/p, we obtain for 
slow rotation: 



AJ = 



544 ,G 



tt'^—pIuR^ = 



34 CM 



(85) 



2835" c^"^" 315 c^R 

From the last two equations we get, after some calculations during which we exploit again 
relationship (68), the following expression for the density functional r[p]: 

"R 



r[p\ 



63 



1 



17poi?2/^«rV(r)(i' 



47rG J, 



uir) + 



P^ 
p{r) 



r p{r)dr 



(86) 



+ 



nR pr ^ pR pr r-r' 

2 / r^p{r)dr I r'^p{r')dr' — / rp{r)dr / r'dr' / r"^p{r")dr" 
Jo Jo Jo Jo Jo 



4.6. The density functionals for the PN transformation functions di and d2 

In eq. (46) we have introduced the two density functionals ii[p\ and when consid- 
ering the transformation from conformal radial coordinates to Newtonian radial coordinates, 
which is ruled by cq. (11) in the case of slow rotation. Therefore the expressions for p[p] 
and uIp] have to be derived from those for di and d2- 

Since we are interested only in PN corrections of order C(c~^), we retain the full ex- 
pression (12) for the transformation function di but only the last term of expression (13), of 
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course evaluated in the case of rigid rotation, for ^2- Therefore on the surface of the spherical 
star we have: 



di{R) = G 
d2{R) - 



m{R) 1 



R 



1 /-^ \ 
+ — j m{r)rdr \ 



(87) 



r dr 



Considering first the function di, exploiting eq. (59) we can write: 



di{R) = 



AttG 



/ r^p(r)dr + — / rdr / r'^p(r')dr' 
Jo R Jo Jo 



(89) 



which for constant density gives, after some calculations: 

d,{R) = -nGpoR' = Ig— 
5 R 

The expression for the density functional ii[p\ thus results, from the above equations: 



(90) 



n 5 1 



2poR^ 



pR pR pr 

/ r^p(r)dr H / rdr / r'^pir')dr' 

Jo R Jo Jo 



(91) 



On the other hand, moving to the transformation function ^2, we can see that eq. (88) 
is independent of the density distribution p(r), and thus it is always: 



This implies the simple identity: 



^ ' 15 



v\p\ = 1 



(92) 
(93) 



4.7. The density functionals for the mass transformation 

From eq. (14) it is possible to obtain the relationship between the baryon rest mass Mq 
and the conformal mass of the fluid configuration. In effect, exploiting eq. (59) we get: 



pR. C ( f^' 

4tt R'^p{R')dR' + i8TT^^[ R'p{R')dR' R"^p{R")dR" 
Jo c \Jq Jo 

+ r R'^p{R')dR' r R"p{R")dR"] +^n^ ^ R'*p{R')dR' (94) 
Jo Jr' J 3 c 7o 
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Introducing then three new density functionals and the conformal mass Mc — (47r/3)po-R^, 
in terms of the angular momentum J we can write: 

M. « MM + f ^Cl*(e.O + l^^^fM (95) 

where the shape functions given in eq. (26) have been added in order to generahze to the fast 
rotation case, as done in §3.3 for the PN total energy, with the same criterium in choosing 
the numerical factor of the PN kinetic correction. 

From what we have learned up to now about the calculation of the explicit expressions 
of the density functionals, it is easy to verify that the two new functionals ([p] and ri[p] 
coincide respectively with the already treated (3[p] and 7[p] (in the case of an application 
of Fubini's theorem for repeated integrals is required). Considering the other new entry, we 
obtain: 

^^P\ = -^^[ ^'P^^)^^ (96) 



5. Analytic determination of the PN onset point of instability 

Now we are ready to extend to PN configurations the treatment of the bar mode in- 
stability made by BR in the Newtonian case and briefly described in §2. We will follow 
strictly their energy variational method, but adding the PN terms in order to determine 
analytically the critical value Cc at which the nonaxisymmetric Jacobi sequence bifurcates 
from the axisymmetric Maclurin one in the case of PN arbitrarily fast rotating stars. 

First, we make for PN conflgurations the same assumptions (i) and (ii) reported and 
discussed in §2. Then we write the total energy of the fluid as in eq. (54): 

E^W + K + U + i^W + AX (97) 

We will consider the total energy, for a given baryon mass Mq, as a function E = E{V, S, J) 
of the conformal volume, the entropy and the angular momentum. The gravitational and 
rotational energies W and K, and the PN gravitational and kinetic corrections AW and 
AK are given in this order from eq. (55) to eq. (58). Their dimensional forms are expressed 
in conformal coordinates, while the shape functions g, /, h, I which they contain are given 
respectively in eqs. (6), (5), (52) and (53). The details on all density functionals involved 
in their definitions are given in §4. Considering finally the internal energy U, since for 
assumption (i) it does not depend on the rotating fiuid shape, in the PN treatment its 
explicit form is not needed. 



-27- 



5.1. PN equilibrium configurations 

In order to construct the equihbrium sequence of relativistic rotating configurations, we 
must minimize the total energy E keeping fixed the baryon mass Mq and not the conformal 
mass M. Therefore we require the validity of both equilibrium conditions: 



dE 


dEdR , 
dR~de ^ ' 


(dE 


de 


\^~de 


dE 


dEdR , 
dR~d^^^ 


(dE 







(98) 



R 



= (99) 

R 

at e 7^ 0, ^ = 0, with the constraint cLMq = 0. This constraint permits us to obtain the 
variation of R, since: 

SM. aM„aR^(9^^ ^^^^ 

R 



de dR de \ de 
gives for the polar eccentricity: 



— = K^- — — 101) 



' ^ 

Recalling now eq. (95), at PN order we obtain that: 

^ _6GMC[p] _ 5 J' vlp] 

and similarly for the variation with respect to the equatorial eccentricity S,- By factoring out 
the parts of the energy contributions independent of e, ^ (and calling them Wq, Kq, AWo, AKq), 
the two equilibrium conditions at PN order give: 

^{Wog + Kof)Re + Wo9e + Kof, + AWoK + AKole = (103) 

rt 

^{Wog + Kof)R^ + Wog^ + Kof^ + AWoh^ + AKok = (104) 

Jrt 

The combination of these last two equations gives, after rearrangement: 



Wo \ge g^J Wo \ge g^J Wo \ge g^ 

(Ko \ 25 _J^v[p] f fe /A . ..... 

^W-VT2^M^^m[vrvj^' ^'''^ 

Now, by recalling definitions (5)- (6) and (52)- (53) respectively for the shape functions 
/, g, h, I, it is possible to verify the identities: 

hm{f,g^-gj^) = (106) 
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lim{heg^ - Qeh^) = 

liMleQi - 9ek) = 



(107) 

(108) 



which solve eq. (105). It must be pointed out that this solution to eq. (105) has an 
important physical meaning: there is always an equilibrium rotating configuration for any 
polar eccentricity e ^ but no equatorial eccentricity (^ = 0), independent of the EOS, 
which in eq. (105) is represented by the three ratios Kq/Wq, ^Wq/Wq and A/Tq/M^o- 

Introducing now in eq. (103) the explicit expressions of the factors independent of e,^ 
(identified by the subscript 0), we obtain: 



5 

'r 



5 J2 
+- 



7[p] „ UG'M' 

\Ie + 



AME?a'^[pY^ 85 c^R^ 



6/?[p]/^[p] + Y4^W ) he 



(109) 



llbGJ'' 1 /85 7[p]r[p] 5^. . ^ . r i r i 23 ^ . , , 



Combining this expression with eq. (102), we get: 



7[p] „ ^ GM 

Je I 



3 5_P 

'5^^^^^'^ AGM^Ra^lpY^ ' c^R 

15 j2 mi[p] 



5 e\p] 



99e + 



5 J' mi[p] 

AGM^R e\p] 



gfe 



+ 



2 GM^Ra^pplp] 



f9e + ^ (gPIpMp] + he 



175 J2 



.seGM^Ra^p^ yi"-^ - ^^[^Hp] - iHpMp] + ) k 



(110) 



+ 



125 



48 Gm^R^a''[p]e[p] 



ffe 



- 



In order to rccxpress this latter equation in terms of the adimcnsional parameter O^/ (nGpo), 
with po mean density of the configuration, we first combine eqq. (40) and (46), thus obtain- 
ing, for slow rotation: 



J 



;MQi?V[p] 



1 + 



GM fl2 



c^R 



- 634? 



(Ill) 



Generalizing to arbitrary fast rotation with SZ's shape functions given in eq. (23), we can 
write: 

35 GM /12 . , 17r[p]' 



1 + 



82 c^R 



63a[p] 



P3(e,0/(e,0 



(112) 
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where the expression of the function ps in terms of polar and equatorial eccentricities is: 

(113) 



65 24^3(1-6 



5/ ' 35 / (1-0^/=^ 
Now from eq. (112) we obtain: 

J2 3 



2^2/3 18 (i_^)i/3 

f— (Al+^2)- ^ 



GM^R nGpo 25 f 



35 



(1 - e2)2/3 140 



(1 - e2)2/3 



1 + 



35 GM 
41"^ 



12 r 1 17 ^[P] , , 



(114) 



and introducing this result in the equilibrium condition (110) we finally get the param- 
eter f2^/(7rGpo) as a function of the compactness parameter GM/(c^R) along ellipsoidal 
equilibrium configurations with polar eccentricity e for any matter distribution. We fix at 
GM/{c^R) — 0.150 the end of validity of our PN approximation. Exploiting the fact that 
a\p\ — 1 independently of the EOS (see §4.4), at the PN order this function is: 



14^ 

H he 

85 



^ ^ Afg^m _ 20 ff3[p]GM 
TrGpo fe l[p] 3 fe 7[p] c^R 
1 S\p\ \ 21 gj 



U(3[p]J 134 /e 



21 m, 21 fglm 

5^^'e[p] ^ 5 /e e[p] 

23 a[p] 



85 , , 5 (3[pHp] , , r 1 
— T p - - , , - 15u p + 
63 2 7[p] ^^^J 14 7[p] 



(115) 



21 ^ 



yMp] - ^r[p] 



Its graphical representation in the case of homogeneous ellipsoids {i.e., when all the density 
functional take the value 1) is reported in Fig. 1, where we can see that in the relativistic case 
the value of Q^/ (nGpo) is larger than the Newtonian one at any given polar eccentricity. This 
confirms the results already found by SZ (but see discussion in §7) and by Chandrasekhar 
(1965b). 



5.2. The PN secular instability point 

In order to determine the point of secular instability in the PN case, we must minimize 
the total energy E with respect to the volume V, again keeping fixed the baryonic mass Mq 
and not the conformal mass M. We start introducing the function of S, Mq, J, independent 
of e: 

U = H-U = PV (116) 
The minimization of the total energy E with respect to the volume V gives: 

dE dE dM fdE\ , 



0.2 0.4 0.6 0.8 1 

e 



Fig. 1. — Adimensional ratio /{nGpo) as a function of polar eccentricity c for PN equi- 
librium sequences of homogeneous ellipsoids. The different curves correspond to 6 equally 
spaced values of the compactness parameter GM/(c^R) increasing from to 0.150. The 
Newtonian Maclaurin sequence (dashed line) is that with GM/(c^R)=0. 

while the constraint dMo = imphes: 

M 



dV dM dV V dV 
The variation of M with respect to V is thus given by 

dM _ 

which, recaUing eq. (95), at PN order results: 



(119) 



aM^i^ecMHIPl (1.0) 



dV V \h c^R e[pr 6c'^MR^ a'^[p]e[p]' 

Therefore from eq. (117) we obtain the scalar virial equation for PN configurations in the 
following form: 



3' ' ' '\5c''Re[p]'' 6cm^R^a^[p]e[p] 
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If we consider the asymmetrical case, where ^ 7^ 0, the total energy is a function 
E — E{S, Mq, J,^) which can be expanded in the neighborhood of ^ = and keeping 
constant S, Mq, J, obtaining, similarly to eq. (1): 



E = Eo{S, Mo, J) + ^Ei{S, Mo, J) + ^M^, Mq, J) + ^^E^iS, Mo, J) + ... 



(122) 



Since at constant S, Mq, J we have dIl\s,Mo,j — 0, the virial equation implies: 



d{W + 2K + 2 AW + 3AK + {K - 2W) 

v[p] 



ISGMCM 5 



S,Mo,J 



and furthermore: 
dE\s,Mo,j 



1 1 

20^W^|5,Mo,J - -dAK\s,Mo,J 



v[p\ 



f 



(123) 



(124) 



S,Mo,J 



We now introduce E, the derivative of e with respect to ^ at constant S, Mq, J. From eq. 
(124) we obtain the following expressions for the expansion terms Ei and E2, remembering 
to consider constant the variables that appear as subscripts: 
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(125) 



(126) 
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d_ 
d_ 
d_ 



(dAK\ 



fdAK\ 



R 



v[p] 



oe 



J' V[P\ 



with the derivatives calculated at ^ = 0. 



J R 



R 



To determine analytically the critical eccentricity Cc we must investigate the condition 
E2 — 0. For this, we exphcit in the right-hand side of eqs. (125)-(126) the shape functions 
f,g,L By factoring out the parts independent of e,^, and using again the notation fx — 
df /dx, we can write the PN expansion coefficients in the form: 



E2 1. ^9 ^ ^ _ , / 42 GMClp] . v\p] 

-^(/^^ + 2S/e5 + SVee + /e^^ + ^feK) (128) 

/18^ GMC[p] r^[p] 85 v[p] \ 

\5 Wo^ c^R e[p\ Wocm^R^ (T^[p]0[p] I2^cmm^ (^^[p]0[p] J 

f9GMC[p] 5 J' 4p] 



Now condition (123) allows the determination of E. In fact, separating the shape func- 
tions from the parts independent of the two eccentricities e, ^, after some calculations we 
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obtain: 



^ R^g^^ Wog^ R Wo'' ^ Wo g^^ Wo g^ 
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(130) 
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Now, from the identities (106)-(108) we get: 
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9e 




le 




fe 



4e 



Inserting these values into eq. (127) we can verify that also in the PN treatment it is Ei 
If we define: 

X{t e) = -k/le 
m^) = -h/fe 

from eqs. (129)-(131) we get, at ^ = 0: 
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Eq. (128) can thus be rewritten as: 
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As a consequence of the definitions of 0, % and ^Z', the quantities in braces vanish identically, 
so that in the end we get: 



^2 1 rv. i 1 A 42 GMC\p\ 



J2 



18 7^0 GMC,{p\ 



(141) 



X 



18^0 GMC\p\^K^. J2 
9 /.r i +577r/- 



85 



J2 



5 W^o c^i? ^[p] W^o cm-^B? a\p\Q\p\ VT c^W^B? a'^\p\Q{p\ 
Inserting in the square brackets the corresponding expressions, and making use of the result: 
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(142) 



which derives from equilibrium condition (103) and exploits the fact that a\p\ = 1 indepen- 
dent of the EOS, after some rather lengthy calculations wc finally obtain the full form of the 
PN expansion coefficient E2 in terms of first and second derivatives of the shape functions 
f,9,h,l- 
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(143) 
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The full expressions of these derivative functions are reported in Appendix B. 
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6. Evaluation of the PN critical point 



6.1. 



The Newtonian limit 



We are now able, by using eq. (143), to evaluate the critical eccentricity Cc, where the 
Jacobi nonaxisymmetric sequence bifurcates from the Maclurin axisymmetric sequence, for 
any PN configuration. But first, in order to check that our PN treatment is consistent with 
the Newtonian results obtained by BR, we analyze the simplified expression of £'2/1^0 after 
rejecting all the PN terms. We thus obtain: 



At first glance, both eqs. (144) and (143) are characterized by factors depending upon the 
mass distribution, that is, the ratios Kq/Wq, AWq/Wq, AKq/Wq, and factors depending 
on the eUipsoid shape. But the equihbrium condition (dE/de) — 0, if considered in the 
Newtonian case, imphes that the ratio Kq/Wq which appears in the Newtonian result for E2 
depends only upon the eccentricity e: 



Therefore eq. (144) depends only on shape functions, and the condition E2 = has an 
unique universal solution Cc for any internal structure of the rotating configuration, which is 
Be = 0.81267, as found in BR. 

On the other hand, the addition of the PN terms in the equilibrium condition (103) 
makes a separation between physical and geometrical quantities in eq. (142) impossible. 
Thus the PN result (143) for E2 remains a mixed expression of factors depending upon the 
internal structure and factors which are functions of the configuration shape. Therefore, in 
the case of PN rotating ellipsoids, the critical value Cc which solves the equation E2 — is 
not universal: to obtain this value we must know the physical quantities M, V, J, together 
with the distribution of mass expressed by the density functionals. 




ge9^ 

fe9e 




(144) 



6.2. The case of incompressible fluids 



Moving now to PN configurations, we see from eq. (143) that the dependence of Cc 
on the rotating eUipsoid internal structure is via the three ratios Kq/Wq, AWq/Wq and 
AKq/Wq. RecaUing expressions (55)- (58) and the definition of 0-subscript quantities as the 
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parts independent of the two eccentricities e, ^, we obtain for these ratios: 
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(146) 
(147) 
(148) 



where we have substituted the volume V with the conformal radius R — (3y/47r)^/^. 



By exploiting eq. (142) is then possible to obtain the expression of for any PN 
rotating configuration at equilibrium: 
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(149) 
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By inserting this expression in eqs. (146) and (148), after some calculations in which we 
exploit the facts that a[p\ — 1, 77 [p] = 7[p] and = for any EOS, we can obtain also 
the PN-approximated ratios Kq/Wq and AKq/Wq in terms of the configuration compactness 
parameter GM/ic^R): 
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(151) 



In this way, we see that with eq. (143) we can evaluate the critical value Cc in the PN approx- 
imation for any given mass distribution and any compactness parameter of the configuration. 



We start by considering the case of constant density mass distribution, i.e., the case 
analyzed by SZ. For homogeneous configurations, we already know that all the density 
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functionals discussed in §4 simply take the constant value 1. Considering different values 
for the compactness parameter GM/{c'^R), the equation £^2 = gives different values of 
the critical eccentricity Cc at the secular instability onset point, as reported in the first two 
columns of Table 1 and in Fig. 2. The corresponding critical values for the adimensional 
ratio fl'^/{7TGpo) are given in the first two columns of Table 2. 
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Fig. 2. — Critical eccentricity ec of the PN secular instability point as a function of the 
compactness parameters GM/(c^R), in the case of constant mass density distribution. 

Another indicator for the onset of instabihty is the ratio K/ | 14^ | of the kinetic energy 
to the absolute value of the gravitational potential energy. A relativistic analog of such a 
ratio can be defined as: 



K _ inj 

JW] ~ \^J-E ^^^^^ 

where the total energy E must be expressed as a function of the angular velocity Q by 
exploiting eqs. (41) and (46), and following the same procedure of §3.3. Therefore in 
axisymmetry this parameter, which is gauge invariant for rigidly rotating objects, at the PN 
order results: 
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In this expression they also appear SZ's shape functions pi2 and p^. In terms of polar and 
equatorial eccentricities, the latter is already given by eq. (113), while the former is: 



/ .^ _ 27(7 171 As jl-e^y/^ 111 (1^)V3 

where again the integral is that reported in Appendix C of SZ. The critical values of the 
ratio K/ | | at the onset point of instability for constant density mass distributions are 
given in the first two columns of Table 3. 
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6.3. Compressible fluids: the polytropic case 

In order to evaluate the critical eccentricity Cc for a more general mass distribution, we 
consider now polytropic distributions. 

By exploiting the properties of such equilibrium configurations, which are exhaustively 
described, e.g., in Chandrasekhar (1957), we are able to calculate the values of all the density 
functionals treated in §4 for any polytropic mass distribution, i.e., for any polytropic index 
n, where P = kp^^^^^. To do this, we must consider each density functional, rewrite its 
expressions in the case of a polytropic mass distribution, and then evaluate such expression 
for any index n. We will thus obtain the considered density functional as a function of the 
polytropic index. 

However in the Introduction we pointed out the result of Bonazzola, Frieben & Gour- 
goulhon (1996) concerning the critical polytropic index for the onset of the viscosity-driven 
bar mode instability, reporting that they find a critical value, for very relativistic objects, 
slightly lower than the Newtonian one: n ~ 0.71 versus n — 0.808. This means that for 
intermediate PN configurations the maximum polytropic index for the onset of the viscosity- 
driven bar mode instability lies between these two values. 

Moreover we must point out here that our PN energy variational method does not 
provide by itself a critical value for the polytropic index. In fact, no energy variational 
method can be sensitive to the mass-shedding limit mentioned in the Introduction, which is 
a dynamical phenomenon. 

The papers which have found critical values of the polytropic index n for the onset of the 
bar mode instability are based on dynamical treatments of the rotating configurations, like 
those used by James (1964) and Bonazzola, Frieben & Gourgoulhon (1996). On the other 
hand, the discussion of uniformly rotating equilibrium models beyond the mass shedding 
limit is motivated in Lai, Rasio & Shapiro (1993) with the fact that they are reasonable 
approximations for the interior of the more realistic, differentially rotating structures, which 
can probably exist beyond this limit. We also remark that for n = 0.5—1.0 one obtains models 
with bulk properties that are comparable to those of observed neutron stars (Stergioulas 
1998). 

Because of all the above arguments, we will assume that the maximum polytropic index 
for the onset of the viscosity-driven bar mode instability in strictly rigidly rotating configu- 
rations, is n = 0.8. 

The determination, for each density functional, of the function which gives its values in 
the polytropic index range n = — 0.8, leads to the following expressions (we omit those 
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density functionals whose unitary value is independent of the EOS) : 
3 
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The value of these density functionals for different polytropic indexes n can be calculated 
by numerical integrations of these equations. In Figure 3 the six functionals 7[p], 5[p\, 
a[p], t[p], ii[p] are reported as functions of the index n. 

In the case of we can compare our result with a formula for the potential energy 
of polytropic spherical distributions due to Betti and Ritter (Chandrasekhar 1957, Chapter 
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IV, eq. (90)): 



W 



3 GM^ 



5-n R 

This latter equation, together with our eq. (55), imphes that it must be j3[p] 
which is exactly the function that we find. 



(162) 
5/(5 -n). 



1.4 - 



1.4 



0.5 




Fig. 3. — Density functionals f3[p],^[p],6[p],a[p],r[p], fi[p] as functions of the polytropic index 
n. Note the change in scale of the ordinate axis in the plots for 5[p] and r[p]. 

Now we have all the pieces to evaluate the critical eccentricity Sc where the Jacobi non- 
axisymmetric sequence bifurcates from the Maclaurin axisymmetric sequence, marking the 
PN onset point of the secular bar mode instability, for different polytropic mass distribu- 
tions and for different compactness parameters. Our results are presented in Table 1. For a 
reproduction in terms of the parameters O^/ (ttG^o) and K/ | | , we refer respectively to 
Tables 2 and 3. 

Of course, the density functional introduced for our PN treatment of bar mode insta- 
bility and whose expressions have been determined for the first time in the literature in §4, 
can be evaluated for any compressible fluid, not only for those which can be described by a 
polytropic EOS. 
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Table 1: Critical eccentricity ec at the PN secular instability point for different polytropic 
indexes n and different compactness parameters GM/ (c^R). The column for n=0 corresponds 
to the case of constant mass density distribution. We fix at GM/(c^ R)=0.150 the end of 
validity of our PN approximation. 
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Table 2: Critical value of the parameter Q^/ (nGpo) at the PN secular instability point for 
different polytropic indexes n and different compactness parameters GM/(c^R). The column 
for n=0 corresponds to the case of constant mass density distribution. We fix at GM/(c^ 
R)=0.150 the end of validity of our PN approximation. 
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Table 3: Critical value of the ratio K/|W| at the PN secular instability point for different 
polytropic indexes n and different compactness parameters GM/(c^R). The column for n=0 
corresponds to the case of constant mass density distribution. We fix at GM/(c^ R)=0.150 
the end of validity of our PN approximation. 
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Therefore in this paper we deposit all the instruments needed for the evaluation of 
the critical eccentricity Cc where the nonaxisymmetric Jacobi sequence bifurcates from the 
Maclurin axisymmetric sequence, once given the fluid configuration EOS. 

More realistic equations of state have been considered in the literature, both in the form 
of schematic analytical models and in the form of numerical tables interpolated by means of 
particular functions. By inserting in the general expressions given in this paper the equations 
relative to a specific EOS, or their numerical representations, the onset point of bar mode 
instability for that specific EOS can be evaluated. This may be the object of future work, 
since we know that new analytical functions are in preparation (P. Haensel 2002, private 
communication) to describe the internal structure of real neutron stars. 



7. Discussion 

We can now discuss the results obtained in last section, comparing them, when possible, 
with other works in the literature. 

First, in the uniform density case, treated in §6.2, we can compare our results with those 
obtained by SZ, which are given in Table 3 of their work. The comparison of the results 
is not straightforward, since they adopt a different compactness parameter GM/{c^Rs), 
where M is the total mass-energy of the configuration while Rs is the equatorial radius 
in Schwarzschild coordinates of the spherical, nonrotating configuration. Considering the 
maximum value for their compactness parameter, SZ point out that it is determined by a 
limit imposed by the relationship between the conformal and the Schwarzschild compactness 
parameters in the spherical limit. This function (reported in their eq. (149)) has a maximum 
for GM/{c^Rs) ^ 0.28, corresponding to GMJic^ai) 0.134, and therefore SZ state that 
their PN formalism can be used to investigate relativistic sequences up to a maximum value 
[GM/ {c^ Rs)]max ~ 0.28. Note that the maximum value which we have fixedioi our conformal 
compactness parameter, le., GM/ (c^R) = 0.150, is similar to the value that SZ's conformal 
parameter takes in correspondence with [GM/{c^Rs)]max- 

However, if the range of compactness parameters considered is almost the same, the 
range of critical eccentricities which we present in this work is significantly smaller than that 
found by SZ. Therefore we can state that also with our treatment we find that the critical 
value of the eccentricity for the onset of the bar mode instability increases as the star becomes 
more relativistic, in the regime in which the PN approximation is valid, but such increase is 
less marked than in SZ. Thus the presence of a stabilizing effect due to general relativity on 
the Jacobi-like bar mode instability is confirmed but weakened in its significance. 
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Another difference between our PN treatment and SZ's appears when considering the 
equihbrium sequences of the rotating configurations. Comparing our Fig. 1 with SZ's analog 
Fig. 1, it is possible to note that we find a much marked increase of the angular velocity 
in homogeneous ellipsoids with the compactness parameter, at any given value of polar 
eccentricity. This discrepancy is solved if wc adopt in our calculations the overall numerical 
factors given by SZ (remember the dilemma of §3.3). With these values we obtain the PN 
equilibrium sequences reported in Fig. 4, which results identical to SZ's Fig. 1. 

0.6 I ' ' ' \ ' ' ' \ ' ' ' \ ' ' ' \ ' ' ' 1 




0.2 0.4 0.6 0.8 



Fig. 4. — Same as Fig. 1, hut obtained with the overall numerical factors given by SZ for 
the PN kinetic corrections (see discussion in the text). 

If we consider now the results we obtained in the more general compressible polytropic 
case (see §6.3), from Table 1 it is evident that the increase in the critical value of eccentricity 
for the onset of the bar mode instability with the compactness parameter is confirmed at 
any polytropic index value, in the regime in which the PN approximation is valid. Thus the 
presence of a stabilizing effect due to general relativity on the Jacobi-like bar mode instability 
is a property also of softer equations of state with respect to the incompressible case. 

The only exception to this general trend can be found in the last column of Table 3, 
which gives the critical values of the ratio K/ | | for a polytopic mass distribution with 
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n = 0.8. This column suggests that increasing the mass concentration towards the center of 
a rotating configuration, i.e., increasing the value of the polytropic index n, there is a value 
above which the viscosity-driven bar mode instability is strengthened by relativistic effects 
and no more weakened. We find that such a value lies in the range 0.7—0.8. However, as 
pointed out in §6.3, in the same range also lies the maximum polytropic index for the onset 
of the instability. Therefore the column for n = 0.8 in Table 3 may be non-representative 
of the viscosity-driven instability, and in this case the value of the polytropic index where 
we find the inversion in the instability strength may also define the maximum index for the 
onset of the instability. 

The increase of the critical value of eccentricity with the polytropic index n for a given 
compactness parameter GM/{c^R) is almost neghgible, and therefore the onset point of 
secular instability in our PN treatment may be considered independent of the polytropic 
mass distribution. 

As we have reported in the Introduction, the same result was reported, but limited to the 
Newtonian case, in Lai & Shapiro (1995). On the other hand, the numerical investigation of 
the viscosity-driven bar mode instability carried out by Bonazzola, Prieben & Gourgoulhon 
(1996) found a weak dependence of the onset point of instability on the polytropic index n 

in Newtonian configurations (sec in particular their Fig. 3). This discrepancy may be due to 
the "ellipsoidal approximation" adopted in the analytical works on rotating configurations 
such as Lai, Rasio & Shapiro (1993), Lai & Shapiro (1995), BR and of course ours. Numerical 
investigations do not need such an approximation, and therefore they lead to slightly different 
results on this item. 

Finally, Bonazzola, Fricben & Gourgoulhon (1996) were not able to investigate the 
secular instability in incompressible fiuid configurations, due to the problems that their 
numerical code had in treating the strong discontinuity of the density profile at the surface 
of the star, which varies suddenly from its constant value to zero. Therefore up to now no 
results from relativistic numerical investigations are available on the dependence of the onset 
point of instability upon the configuration compactness. A new numerical code, which solves 
the discontinuity problem (Gondek-Rosinska & Gourgoulhon 2002), makes the comparison 
possible between a fully relativistic numerical investigation and SZ's and our analytical PN 
treatments. 
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8. Conclusions 

In this paper wc have treated the bar mode secular instabihty of rigidly rotating equilib- 
rium configurations for neutron stars with an analytic energy variational method in the PN 
approximation. The method, which is the extension to PN configurations of that used by 
BR for the Newtonian treatment of bar mode instability, gives results for any (but defined) 
EOS. 

After the full derivation of the equation which gives as its solution the critical eccen- 
tricity Cc at which the secular bar mode instability sets on, we considered some particular 
equations of state and evaluated the corresponding critical values. To do this, we introduced 
density functionals which allowed the generalization of the physical quantities involved in 
the treatment from the constant mass density to the arbitrary density profile form. The 
determination of the explicit expressions of such functionals has been made for the first time 
in the literature. 

We started by checking that in the Newtonian case, i.e., when all the PN corrections are 
null, the critical value for the eccentricity is Cc = 0.8127, as found in the precedent literature 
on this item. 

Then we considered PN configurations with a constant density mass distribution. In 
this case we have found that the critical value Cc depends on the neutron star compactness 
parameter GM/{c^R), resulting larger for more relativistic, i.e.. more compact, objects. 
Thus the effect of general relativity, considered in the PN approximation, is to weaken the 
bar mode instability, stabilizing the object against secular transition from the axisymmetric 
Maclaurin to the triaxial Jacobi sequence. This is consistent with the results obtained by SZ 
in their work on incompressible rotating stars, but we find a less marked stabilizing effect. 

Finally we considered polytropic mass distributions. For these configurations we have 
calculated numerically the values of all the density functionals involved in our PN energy 
variational method, obtaining them as functions of the polytropic index n up to the dynam- 
ical limit n ~ 0.8 imposed by the mass-shedding of the rotating configuration. The result is 
the confirmation of the increase of the critical eccentricity with the compactness parameter 
also for n > (softer) polytropes, of course in the regime in which the PN approximation is 
valid . 

This latter investigation has also shown that for a fixed compactness parameter the 
increase of the critical value of eccentricity with the polytropic index is negligible, thus 
extending to PN configurations the independence of the onset point of secular instability 
upon the polytropic mass distribution, at least in the "ellipsoidal approximation" regime. 
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The formula for the PN total energy (54), with the expressions of its PN corrections 
(57)-(58), the explicit general forms of the density functionals introduced in our energy 
variational method (62), (69), (77), (78), (82), (86), (91), (93), (96), together with those 
specialized to the polytropic case (155)-(161), and the full expressions of the ellipsoidal 
shape functions together with their derivatives (Appendix B), will all be found useful for 
future investigations. 
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The dimensionless coefficients Ai are given in eqs. (3.33) — (3.35) of Chandrasekhar 
(1969a) in terms of standard incomplete elliptic integrals involving only the values ai{i = 
1,2,3) of the three semiaxes of the ellipsoid outer surface. It is possible to calculate them 
(as in SZ) in terms of the axial ratios Ai = (aa/ai)^/'^ and A2 = (03/02)^^^ and therefore as 
functions of the eccentricities e and ^. The standard incomplete elliptic integrals involved in 
their definition thus become: 
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A. The coefficients A^ in terms of the eccentricities e, ^ 
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After some algebraic manipulation, we obtain the expressions: 

Ai = ^ ; ^ vicO (A5) 

t(e, Oi(e, - 2v{e, Qe^l - Q V^(l - e^) - 2e(l - e^) 

e(e2 - ^) 

where we have defined the auxiliary functions: 

/-arcsine / t \ —1/2 

^(e,^) = — sin^x ( 1 - -^sin^x j (A8) 



/■arcsine / ^ \ ""^^^ 

j(e, — J (l ^sin^xj dx 



(AlO) 



/■axcsin e / t 

j(e,0 = M-^sin^xj dx (All) 

It must be noticed moreover that the coefficients are not independent, since it is valid 
the relationship Ai + A2 + A3 — 2. 



B. The full expressions of the derivative functions in equation (143) 

As we have shown in the paper, the critical value of eccentricity for the onset of bar 
mode instability is given by equating to zero expression (143), evaluated in the limit ^ ^ 0. 
The dependence of this expression on the eccentricity e is contained in the combination of 
derivatives of the shape functions /, g, h, I. We report here the full expressions of these 
derivative functions: 

lim — = lim — = lim — = lim — = (Bl) 

^-^0 Qe 3eVl - e2 + {2e^ - 3) arcsin e 
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9-4e2 27- 30e2 + 4e4 

lim Qef — TT, ^TTT7 77 T^TTiTr aiCSin 6 ( B5 ) 

^-^0^^^ 12e3(l-e2)V3 366^(1 - e2)5/6 ^ ^ 

27+10e2 2 81-24e2-40e\ oxi 
= 96i^(^-')' + 288i^ (l-e)^arcsme (B6) 

lim/.., = 56e9(i ^_ ^2)2/3 [-9e^(-36 + 69e^ - 49e^ + 16e^) 

+2eVl - e2(-324 + 513e^ - 468e^ + 152e^) arcsin e (B7) 
+ (324 - 729e^ + 936e^ - 604e^ + 96e^) (arcsin e)^] 

e^S^^^ ^ ii^^^[9e2(-387 + 615e^-304e^ + 76e«) 

+2eVl-e2(3483 - 43746^ + 3600e^ + 272e^) arcsin e (B8) 
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i^o 140ei3 L ^ ' 
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lim-^ = {9 re2(-225 + 474e^- 289e^ + 34e^-30e^ + 36e^°) 

+2eVl-e2(225 - SQQe^ + 186e^ - 9e^ + 16e^) arcsine 

+9(-25 + 61e^ - 48e^ + 12e^) (arcsin e)^] } (B12) 
/ |70e^(l - e^fl^ [3ev^r^+ (-3 + 2e') arcsine] } 
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